Thermodynamics of strongly interacting fermions in two-dimensional optical lattices 
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We study finite-temperature properties of strongly correlated fermions in two-dimensional optical 
lattices by means of numerical linked cluster expansions, a computational technique that allows one 
to obtain exact results in the thermodynamic limit. We focus our analysis on the strongly interacting 
regime, where the on-site repulsion is of the order of or greater than the band width. We compute 
the equation of state, double occupancy, entropy, uniform susceptibility, and spin correlations for 
temperatures that are similar to or below the ones achieved in current optical lattice experiments. 
We provide a quantitative analysis of adiabatic cooling of trapped fermions in two dimensions, by 
means of both flattening the trapping potential and increasing the interaction strength. 
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Recent optical lattice experiments have opened a new 
venue for exploring the effects of strong correlations in 
quantum lattice models. For example, the superfluid to 
Mott-insulator transition for bosons has been observed 
in geometries of three 

0, two 0, and one [3f dimen- 
sion. Currently, there is a race to access temperatures 
low enough for the transition to the antiferromagncti- 
cally ordered Neel state in three dimensions, or possibly 
more exotic states in two dimensions, to be observed for 
fermions [1,01 ■ So far, the interaction strength and the 
temperature in lattice fermion experiments remain rel- 
atively high in comparison to the hopping amplitude t. 
This is in part because t, which is set by optical lattice 
parameters, is in general small in the regimes where one- 
band models are applicable. 

On the theoretical side, there is an ever-increasing de- 
mand for precise numerical results for the relevant pa- 
rameters of the Hubbard model and for large enough 
system sizes, which could be used to interpret current 
experiments and also provide suggestions for future ex- 
periments For this model, especially for strong 
interactions, the present computations become particu- 
larly challenging as the temperature is lowered below the 
hopping amplitude. 

Here, we study various thermodynamic quantities such 
as the equation of state, entropy, double occupancy, and 
spin correlations in the thermodynamic limit for inter- 
actions up to three times the band width, utilizing nu- 
merical linked cluster expansions (NLCEs) Wc 
obtain a detailed understanding of the evolution of var- 
ious quantities with adiabatically increasing interaction 
strength, of great interest to current optical lattice exper- 
iments. Using the local density approximation (LDA), we 
analyze the thermodynamics of fermions in a harmonic 
trap and calculate their temperature as a function of the 
interaction strength and total entropy. We also present 
a quantitative analysis of various cooling schemes for the 
experiments [15l - fl7l | . 



We consider the two-dimensional (2D) Hubbard Hamil- 
tonian, 



^ = -<E(i^+H.c.)i[/^ 
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(1) 

where c\ a creates (annihilates) a fermion with spin 
a on site i, and hi a — & ia Cia is the number operator. (..) 
denotes nearest neighbors (NNs), U is the strength of the 
on-site repulsive interaction, and Vi is a space-dependent 
local chemical potential, t = 1 (H = 1 and ks = 1) sets 
the energy scale throughout this paper. 



III. COMPUTATIONAL APPROACH 

In linked-cluster expansions [l2| . we express an exten- 
sive property of the model per lattice site in the thermo- 
dynamic limit (P) in terms of contributions from all the 
clusters, up to a certain size, that can be embedded in 
the infinite lattice: 



P 



L(c)w p (c), 



(2) 



where c represents the clusters. This contribution is pro- 
portional to the weight of each cluster for that property 
[iUp(c)] and to its multiplicity [£(c)]. The latter is defined 
as the number of ways in which that particular cluster 
can be embedded in the infinite lattice, per site. The 
weight, on the other hand, is calculated recursively as 
the property for each cluster \P{c)} minus the weights of 
all its subclusters: 



,(c)=P(c)-£> p (s). 



(3) 



Here, we use the NLCE, where V{c) is computed by 
means of full exact diagonalization fl3| . Because of 
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the exact treatment of individual clusters in the NLCE, 
the series converge at significantly lower temperatures 
in comparison to high-temperature expansions in which 
perturbation theory is used [l3j . 

NLCEs are complementary to quantum Monte Carlo 
(QMC) approaches, such as the determinantal QMC 
(DQMC) [IH, or dynamical mean-field theory [l| and its 
cluster extensions, such as the dynamical cluster approx- 
imation (DCA) [H HJ. They can also help to bench- 
mark future experiments as well as new computational 
techniques. This is because NLCEs do not suffer from 
statistical or systematic errors, such as finite-size effects, 
and, as opposed to the DQMC and DCA, they are not 
restricted to small or intermediate interaction strengths. 
In Ref. [IH, we make our raw NLCE data for a wide 
range of interactions available for comparison. 

The validity of NLCEs, however, is limited to a region 
in temperature in which the series converge (the con- 
vergence region). We have found that, for the Hubbard 
model, NLCEs converge down to lower temperatures as 
the strength of the interaction is increased. At half- 
filling, and for interactions larger than the band width, 
NLCEs can access the region with strong antiferromag- 
netic (AF) correlations, identified by the suppression of 
the uniform susceptibility. Although, the method does 
not have any systematic restriction away from half-filling, 
in the latter region, the series fail to converge at temper- 
atures as low as those accessible to the half-filled case. 
This prevents us from accessing low-temperature phases, 
such as d-wave superconductivity, that arguably exist in 
this model at finite doping. 

We begin our analysis with the homogeneous system 
(Vi = 0) in the grand canonical ensemble. For each U, 
we compute all properties for a very dense grid of chem- 
ical potential (//) and temperature, so that we can also 
follow properties at constant density (n) Q ■ The NLCE 
calculations are carried out on the square lattice up to 
the ninth order in the site expansion (nine sites) . We use 
Wynn and Eulcr algorithms for summing the terms in the 
series to extend the region of convergence (l3j . Since only 
NN hopping is considered, all properties of the particle- 
doped system can be expressed in terms of those from 
the hole-doped system. Hence, away from half-filling, we 
only show results for the hole-doped system. 



IV. RESULTS 
A. Equation of State 

The equation of state for the Hubbard model provides 
important information about correlation effects as the 
strength of the on-site interaction is increased, and can 
be studied in optical lattice experiments. In Figs. QJa)- 
(c), we depict the equation of state at three different 
temperatures, T = 0.82, 0.55, and 0.25, for the weak- 
, intermediate-, and strong-coupling regimes (U = 4, 
8, and 12, respectively). For the last two values of U 
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FIG. 1. (Color online) Top: Equation of state for (a) [7 = 4, 
(b) U = 8, and (c) U — 12 and at three different temperatures. 
Except for U = 4 at T = 0.25, NLCE results converge for all 
the values of chemical potential presented here. Only the 
last order of the series is shown after using Wynn sums with 
three cycles of improvement. Bottom: Normalized double 
occupancy vs T at four hole dopings for (d) U = 8 and (e) 
U = 16. We use Euler sums for the last six terms at half-filling 
and Wynn sums for tl/ 1. Thin (black) lines in (d) and (e) 
are the results for the one to last order of NLCEs after the 
above sums. The inset in (e) magnifies the low-temperature 
region for U = 16. The unit of energy is set to the hopping 
amplitude t. 



[Figs. [Tf b) and QJc)], one can see the emergence of an 
incompressible region around ji = U/2, a clear signature 
of the Mott gap opening in the density of states at low 
temperatures. 



B. Double Occupancy 

In Figs. [Ud) and Die), we show the double occu- 
pancy, D = (fifh^,), normalized by its uncorrelated high- 
temperature value (n 2 /4) for U = 8 and 16, respectively. 
The double occupancy exhibits a clear low-T rise with 
decreasing temperature. This feature has attracted a 
lot of attention recently, especially after the real-space 
DMFT study of the three-dimensional (3D) version of 
the model in a harmonic trap |{|. Gorelik et al. argued 
that the onset of the AF ordering in the strong-coupling 
regime is signaled by an enhanced double occupancy, 
which can be directly measured in optical lattice exper- 
iments. However, according to Figs, [lid) and [He), the 
low-temperature rise occurs not only at half-filling, but 
also away from it. Moreover, the rise starts at even higher 
temperatures for higher dopings. This implies that the 
enhancement of D in the trap upon lowering the temper- 
ature is significant in the Mott-insulating core as well as 
in other areas of the trap where density is < 1. There- 
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FIG. 2. (Color online) Entropy (a) vs T at half-filling, and (b) 
vs n at T = 0.82 for different values of 17. Thick (thin) lines 
are results for the last (one to last) order of the expansions as 
explained in Fig. [T] In (a) , we have also included the entropy 
for the AFHM with the exchange interaction J = 0.25. 

fore, in real experiments, such an enhancement can be 
observed in systems that have a very small or even no 
Mott insulating region at the center of the trap. Hence, 
the observation of an increase in D alone may not sig- 
nal the onset of AF order. To ensure that AF order is 
emerging, one must also make sure that the density is 1 
in most of the trap. 

For large values of U [see, e.g., U = 16 in Fig. [He)], 
the normalized D is almost independent of doping below 
T ~ 1 and down to the lowest accessible temperatures for 
n > 0.85 (see inset), implying that Dan 2 in this region. 
One can understand the latter from the fact that local 
moments are likely ordered, and the double occupancy 
arises from virtual hoppings to NN sites, so a relatively 
small number of extra holes only modifies the probability 
of those hoppings (accounted for by n 2 ), not the actual 
process. 

C. Entropy 

Generally, when using QMC-based methods, entropy 
calculations involve numerical derivatives and/or integra- 
tion by parts [HI, 01 , which can introduce systematic er- 
rors. Within NLCEs, the entropy is computed directly 
from its definition in the grand canonical ensemble: 

S = HZ) + ^1^, (4) 

where Z is the partition function. 

We first study the entropy at half-filling. Results are 
shown in Fig. Eta) as a function of the temperature for 
U = 4, 8, and 16. There are two distinct regions of fast 
decrease in the entropy in the strong-coupling regime, 
e.g., U — 16. Those regions are separated by a crossing 
point of curves for different values of U around T = 0.6, 
corresponding roughly to S = ln(2). The emergence 
of these two regions results from the fact that as U 
increases, charge degrees of freedom are suppressed at 
higher temperatures due to the higher price of double oc- 
cupancy, and at the same time, the characteristic energy 




FIG. 3. (Color online) Interaction dependence of different 
quantities at half-filling. Left: (a) Entropy and (c) NN spin 
correlations at fixed temperatures. Right: (b) Temperature 
and (d) NN spin correlations at constant entropies. In (b), 
T* represents a crossover temperature to the region where AF 
correlations grow exponentially with decreasing temperature 
(shaded area). 



scale of the spin degrees of freedom, J = At 2 /U, becomes 
smaller, pushing the low-T drop to lower temperatures. 
We find that in the latter region, the entropy curves for 
large U (> 14) follow very closely the entropy of the anti- 
ferromagnetic Hciscnbcrg model (AFHM). This is shown 
for U — 16 in Fig. [U[a), where we also plot the entropy 
of the AFHM with J = 0.25. 

In Fig.^b), we show the entropy away from half-filling 
for a range of interactions at a fixed T = 0.82. In the 
weak-coupling regime (e.g., U = 4), the entropy increases 
monotonically with the density and is maximal at half- 
filling. Since correlations play a small role, the system be- 
haves similarly to a noninteracting system; i.e., the closer 
to the point where there is an equal number of electrons 
and holes, the higher the entropy. This trend changes 
upon increasing U , for which the moment ordering sup- 
presses the entropy significantly close to half-filling. As 
a result, there is a maximum in the entropy in the vicin- 
ity of ri ~ 0.85 for all interactions in the strong-coupling 
regime [25[. Below, we discuss how these features are 
reflected in the properties of trapped systems. 

We further take advantage of the fact that, within NL- 
CEs, arbitrary values of U can be studied at no additional 
computational cost and determine the dependence of the 
quantities of interest on U. In Fig. EJa), we show S at 
half-filling as a function of U at fixed temperatures. The 
temperature regions identified for the entropy in Fig.^a) 
are more clearly seen in Fig. [3ja); i.e., the entropy stays 
more or less the same for T ~ 0.6 (the crossing point) 
as U increases, whereas it generally decreases (increases) 
with U for T > 0.6 (T < 0.6). 

The possibility of adiabatic cooling with increasing U in 
optical lattices has been studied by a number of groups 
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FIG. 4. (Color online) (a) Uniform spin susceptibility and (b) 
NN spin correlations at half-filling vs temperature for different 
interactions, x peaks at T* , below which AF correlations grow 
exponentially with decreasing temperature. S zz also shows a 
sharp increase around T* . 



recently @, HH, [26|. Considering interactions that arc 
often no larger than 3/2 of the band width, they argue 
that this process is rather weak in two dimensions. Here, 
we revisit the problem by including results for values of 
U up to three times the band width [see Fig. EJb)] . We 
find that if we start at relatively high values of entropy 
below ln(2) (e.g., S = 0.68), accessible to current 2D 
experiments, and continue to increase U, the half- filled 
system can be cooled down to very low temperatures. 
For lower entropy, e.g., S — 0.4, one can even access 
the region with exponentially large AF correlations below 
the crossover temperature, T* , as shown in Fig. Elb). 
We take T* as the temperature where the uniform spin 
susceptibility (x) as a function of temperature peaks [2j . 
We find that T* also coincides with the temperature at 
which the NN spin correlations, S zz = \{J2uj) S Z S Z )\, 
show rapid growth with decreasing temperature. The 
uniform spin susceptibility and the NN spin correlations 
are depicted in Fig. |4j For U » 1, T* is expected to 
scale with the AF exchange constant (J) in the effective 
Heisenberg model, i.e., oc 1/U. This behavior, which 
is demonstrated in Fig. [3] (b), is observed not only for 
T* , but also for the large- U tails of the isentropic curves 
when S < ln(2), as J is the only energy scale in that 
region. For U < 12 in the weak-coupling regime, T* is 
not accessible to our NLCE. Therefore, we have taken 
T* from the DQMC results in Rcf. for that region. 
It is interesting to see that this crossover temperature 
peaks when the value of the interaction is around the 
band width. 



D. Nearest- neighbor Spin Correlations 

What is perhaps more important from the experimen- 
tal point of view is how AF correlations change during 
the process of adiabatically increasing U. As mentioned 
in Sec. [J one of the current main goals in cold fermion ex- 
periments is to achieve AF in the Mott-insulating state. 
However, the challenge in this case lies not only in real- 
izing such a state but also in detecting it. Very recently, 
experimental breakthroughs have been reported which 
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FIG. 5. (Color online) (a) Density, (b) entropy, (c) NN spin 
correlations, and (d) double-occupancy profiles of fermions in 
a harmonic trap with p — 22.9, governed by the Hubbard 
model with U — 16 at T = 0.76. The average entropy per 
particle is s = 0.56. £ = (2dt/V) 1 ^ 2 is the characteristic 
length. 



allow the detection of NN spin correlations [221, [28| . 

NN spin correlations, S zz , can also be computed ex- 
actly using NLCEs. As expected, we find that S zz is 
largest at half-filling for all interactions. Therefore, we 
focus on the half-filled system and plot this quantity per 
site vs U at constant temperatures in Fig. (3^c) and at 
constant entropies in Fig. [3jd). The dependence of S zz 
on U at constant T is nontrivial. As the temperature is 
lowered to T ~ 0.3, a peak develops in the spin correla- 
tions around {7 = 8, which is indicative of the largest ef- 
fective exchange interaction between NN spins. The peak 
is a result of the interplay between weak moment forma- 
tion in the weak-coupling regime (U < 8) and the 1/U 
decrease in the effective J in the strong-coupling regime. 
We find that at lower temperatures (T — 0.21), the max- 
imum of S zz occurs at U ~ 9, which is not expected to 
change significantly with further decreasing temperature. 

At constant entropy, on the other hand, this picture 
is strongly modified. Figure [3Jd) shows that S zz satu- 
rates to a finite entropy-dependent value with increasing 
U along the isentropic paths in Fig. (3jb), provided that 
S < ln(2). Note that, even though adiabatic cooling may 
not be efficient to arrive at regions with large AF corre- 
lations in 2D [23|, the value of the NN spin correlations 
will be maximal in the large-{7(^ 12) region if S < 0.6. 
This is convenient for experiments in optical lattices for 
which U is typically large compared to the band width. 



E. Trapped Systems 

To make direct contact with experiments in optical lat- 
tices, we study the manifestation of our previous results 
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FIG. 6. (Color online) Same as Fig. but for p = 10.8. The 
average entropy per particle is s = 0.85 in this case. 



in systems confined by a spatially varying harmonic po- 
tential, Vi = Vrf. Here, denotes the radial distance 
of each site to the center of the trap, and for any given 
value of U, all properties of the system are determined by 
the characteristic density p = N(V/2dt) d / 2 [2{|, where 
d is the dimensionality and TV is the number of parti- 
cles. The resulting inhomogeneous Hubbard model is 
then studied using the LDA along with our results for 
the infinite system. A recent QMC study of the inhomo- 
geneous Hubbard model [ll[ has shown that the LDA is 
a good approximation for local observables at the tem- 
peratures accessible here. We should stress that NLCEs 
are ideal for this kind of study because, for each value of 
U, one can compute all properties for a very dense grid 
of temperatures and chemical potentials at almost no ad- 
ditional computational cost. The same is, of course, not 
true for QMC-based calculations, where each tempera- 
ture and chemical potential requires a separate compu- 
tation. 

In Fig. [5{a), we plot the resulting density profile for 
U = 16 at T = 0.76. We have chosen p = 22.9 such that 
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FIG. 7. (Color online) Temperature vs characteristic density 
at constant entropies per particle for (a) U = 8 and (b) U = 
16. 




U U 

FIG. 8. (Color online) Temperature vs U at constant en- 
tropies per particle for (a) p — 10.8 and (b) p — 22.9. The 
shaded area, the same as in Fig.[3jb), is the region of exponen- 
tially large AF correlations below T* in the Mott insulating 
core of the trap. 

there are band-insulating (n = 2) and Mott-insulating 
(n = 1) domains in the trap. Very useful information for 
the experiments is provided by the spacial distribution 
of the density, entropy, NN spin correlations, and double 
occupancy, as shown in Fig. [S] The entropy is minimal 
(0) in the band insulator, peaks at n ~ 1.18 and 0.82, 
consistent with Fig. H^b), and has a local minimum in 
the Mott ring. In the latter region, spin correlations are 
maximal, and as expected for this large value of U , the 
double occupancy is large only in the region where n > 1. 

In Fig. [SJ we show the same quantities as in Fig. [SJ for 
the reduced p of 10.8 at the same temperature and inter- 
action strength. As a result of this isothermic change, the 
entropy per particle increases from 0.56 to 0.85. Never- 
theless, the Mott-insulating region with a relatively uni- 
form entropy profile is clearly seen over most of the trap. 

It has become apparent in experiments with fermions in 
optical lattices that cooling approaches beyond the stan- 
dard evaporative cooling techniques are required if one 
is to reach temperatures low enough that exotic physics 
emerges. Three recent proposals have shown how to gen- 
erate low-entropy states where a large fraction of the 
system is in a band- insulating domain (i.e., with large 
values of p) [l6|, [It], H(| ■ The idea is then that one can 
adiabatically reduce the trap strength (the characteris- 
tic density p) so that the effective temperature of the 
fermions decreases. In this way, antiferromagnetism and 
other low-temperature phenomena can be explored. 

In Fig. we show quantitatively how this idea works 
for trapped 2D systems. We plot the temperature as a 
function of p for various values of the total entropy per 
particle. Recent studies have shown that the entropy per 
particle (s) for a particular U can be estimated by fit- 
ting the double-occupancy measurements at different p 
to data from numerical simulations In Figs.JTJa) and 
[7Jb) one can see that, for the two values of U shown, the 
temperature decreases rapidly with decreasing p, demon- 
strating that this approach works very efficiently for 2D 
trapped systems. The inflection point, shown, e.g., for 
s = 0.9 in Fig. [TJb), is the signature of a large Mott re- 
gion forming in the trap (as seen in Fig. [5]) . This occurs 
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provided the entropy is low enough and for a range of 
characteristic densities that depends on U. AF ordering 
in the Mott core emerges at T* , which, for low entropies, 
can be reached before the Mott insulator is destroyed by 
further flattening of the trap. 



It is also interesting to study what happens to the tem- 
perature of a trapped system as one increases the inter- 
action strength at constant entropy. [Results for homo- 
geneous systems at half-filling are presented in Fig.^b).] 
In Fig. |S1 we show isentropic curves in the T — U plane 
for trapped systems at various entropies and for the two 
characteristic densities, p = 10.8 and 22.9, used in Figs. [5] 
and [6] As expected from the results in Fig. the shape 
and location of the isentropic curves depend strongly on 
the value of p. In Fig. [5J we also show the same shaded 
area as in Fig. [3jb) below T*, which, here, represents 
the region where the Mott-insulating core of the trap 
develops large AF correlations. Our calculations show 
that cooling can take place in trapped systems as the 
interaction increases. The entropies at which cooling is 
observed, and the values of U at which cooling occurs, 
depends on the characteristic density in the trap. Hence, 
as reported in Ref. [3l[ for 3D systems, adiabatically in- 
creasing the interaction strength can allow experimental- 
ists to reach the temperatures needed to observe the onset 
of (quasi-)long-range AF correlations in a trapped sys- 
tem. Unfortunately, unlike for the homogeneous system 
at half-filling, our NLCEs do not provide access to the 
temperatures relevant to that region for the 2D trapped 



system. 

V. SUMMARY 

In summary, utilizing NLCEs, which, within the con- 
vergence temperature region are free of statistical and/or 
systematic errors and provide exact results in the thermo- 
dynamic limit, we have calculated thermodynamic prop- 
erties, such as the equation of state, double occupancy, 
entropy, uniform susceptibility, and NN spin correlations, 
of the 2D Hubbard model for a wide range of interaction 
strengths and temperatures. Precise data for the entropy 
on a dense temperature grid allowed us to study temper- 
ature and NN spin correlations, relevant to optical lattice 
experiments, as a function of the entropy. Wc find that 
for any S < ln(2), by adiabatically increasing U to very 
large values, the temperature decreases as 1/U and the 
spin correlations saturate to an entropy-dependent value 
beyond U ~ 12. Using the LDA, we have discussed the 
implications of our results for lattice fermions in the pres- 
ence of a confining harmonic potential. In particular, we 
have shown how cooling can be achieved by reducing the 
confinement strength in a system that starts with a wide 
band-insulating domain in the center of the trap, or by 
adiabatically increasing the interaction strength. 
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